Quantum splines 



Dorje C. Brody 1 , Darryl D. Holm 2 , David M. Meier 2 

1 Mathematical Sciences, Brunei University, Uxbridge UB8 3PH, UK 
2 Department of Mathematics, Imperial College London, London SW7 2AZ, UK 

A quantum spline is a smooth curve parameterised by time in the space of unitary transformations, 
whose associated orbit on the space of pure states traverses a designated set of quantum states at 
designated times, such that the trace norm of the time rate of change of the associated Hamiltonian 
is minimised. The solution to the quantum spline problem is obtained, and is applied in an example 
that illustrates quantum control of coherent states. An efficient numerical scheme for computing 
quantum splines is discussed and implemented in the examples. 
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Controlling the evolution of the unitary transforma- 
tions that generate quantum dynamics is vital in quan- 
tum information processing. There is a substantial lit- 
erature devoted to the investigation of the many aspects 
of quantum control [T]. The objective of quantum con- 
trol is the unitary transformation of one quantum state, 
pure or mixed, into another one, subject to certain crite- 
ria. For example, one may wish to transform a given 
quantum state \ip) into another state \<p) unitarily in 
the shortest possible time, with finite energy resource 
[2]-[4]. When only the initial and final states are in- 
volved, many time-independent Hamiltonians are avail- 
able that achieve the unitary evolution \ip) — >• |</>), and 
we simply need to find one that is optimal. However, 
transforming a given quantum state along a path 
that traverses through a sequence of designated quan- 
tum states — > \<j>\) — » \<j>2) —»•••—>' \4> n ) cannot 
be achieved by a time-independent Hamiltonian. To re- 
alise this chain of transformations in the shortest possible 
time, one chooses the optimal Hamiltonian Hj for each 
interval \<j>j) — » [3JH], and switches the Hamilto- 

nian from Hj to Hj+i when the state has reached 
However, instantaneous switching of the Hamiltonian is 
in general not experimentally feasible. 

In the present paper, we consider the following quan- 
tum control problem: Let a set of quantum states |0i), 
|^2) 9 ■ • • , \<i>m) an d a set of times ti, t<z, • • • , t m be given. 
Starting from an initial state |^o) at time to = 0, find 
a time-dependent Hamiltonian H(t) such that the evo- 
lution path \ip t ) passes arbitrarily close to \<pj) at time 
t = tj for all j = 1, . . . , m, and such that the change in 
the Hamiltonian, in a sense defined below, is minimised. 
The solution to this problem will generate a continuous 
curve in the space of quantum states that interpolates 
through the designated states, just as a spline curve in- 
terpolates through a given set of data points. We thus 
refer to this solution as a quantum spline. 

There is a difference between a classical spline curve 
and a quantum spline. In the classical context the solu- 
tion curve passes through a given set of points, whereas in 
the quantum context, a curve on the space of pure states 




FIG. 1: (colour online) A quantum spline for a two-level sys- 
tem. The lower-left initial state and the targets are repre- 
sented by black dots. The variational formulation of the prob- 
lem requires to minimise a functional that measures both the 
cost related to the change of the Hamiltonian, and the amount 
of mismatch between the trajectory and the target points. 



in itself has no operational meaning. Thus, instead of 
finding a curve in the space of pure states where the des- 
ignated states lie, we must find a time-dependent curve 
in the space of Hamiltonians that in turn will generate 
the curve in the unitary transformation group needed to 
produce an optimal trajectory. In other words, we shall 
seek a curve in the associated Lie algebra, which of course 
is equivalent to the space of Hamiltonians, up to multi- 
plication by i = ^/—T- 

Our approach involves variational calculus in the Lie 
algebra of skew-Hermitian matrices, with constraints 
that take values in the unitary group. In addition, since 
our optimality condition for quantum splines involves the 
time-derivative of \H(t), we shall make use of the tech- 
niques developed recently for the higher-order calculus of 
variations on Lie groups and their algebras [5j [6] . By ex- 
tending these results we are able to: (a) derive the Euler- 
Lagrange equations Q and ([9| below that solve quantum 
spline problems; and (b) devise an efficient discretisation 
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scheme to numerically implement the solution. An ex- 
ample of such a solution for a two-level quantum system 
is sketched in Fig. [I] As an application, we illustrate how 
the results transform a quantum state along a path that 
lies entirely on the coherent-state subspace. 

The optimal curve H (t) that solves the quantum spline 
problem is the minimiser of a 'cost functional' (action) 
consisting of two terms: The first term measures the over- 
all change in the Hamiltonian during the evolution. For 
this purpose we shall consider the trace norm, i.e. for a 
pair of trace-free skew-Hermitian matrices A and B we 
define their inner product by 



(A, B) = -2tr(AB), 



(1) 



where the factor —2 is purely conventional. Thus, if H 
is a time-dependent Hamiltonian and H its time deriva- 
tive, the instantaneous penalty arising from changing the 
Hamiltonian is given by \(iH,\H) = tr(H 2 ). The second 
term penalises the 'mismatch' between the state \ipt-) a t 
time tj and the target state \4>j). For this purpose we 
shall use the standard geodesic distance: 
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for a pair of states and \<j>). Writing U(t) for the 
parametric family of unitary operators generated by H(t) 
so that \iptj) = U(tj)\%jjQ), the mismatch penalty is chosen 
to be ^D 2 (U(tj)i/jQ, (j)j)/(T 2 , where the tolerance a > is 
a tunable parameter so that the penalty is high when a 
is small, and the factor of a half is purely conventional. 

The action, of course, must be minimised subject to 
the constraint that the dynamical evolution of the state is 
unitary. That is, U must satisfy the Schrodinger equation 
U = —\HU, in units H = 1. Therefore, given an initial 
state |^o) at time to = 0, a set of target states |0i), 
• • • , \4>m) at times ti, 

H(0) = Hq, we wish to find the minimiser of 



, t m , and an initial Hamiltonian 



J = j™ (§ (Iff, Iff) + (M.UU- 1 +i#))d£ 
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(3) 



where the minimisation is over curves U(t) G SU(n + 1) 
and iH(t), M(t) G su(n + 1). Additionally, we require 
smoothness of these curves on open intervals (£/, tj+i) for 
j = 0, . . . , m — 1; U(0) = 1; and the continuity of U(t) 
and H(t) is assumed everywhere. The curve M(t) acts as 
a Lagrange multiplier enforcing the kinematic constraint. 

Before we proceed to vary the action J let us com- 
ment on the choice of the initial Hamiltonian Hq. We let 
Hq be such that the trajectory e~ lHot \ipo) corresponds 
to the geodesic curve on the space of pure states joining 
|^o) and |0i); the construction of such a Hamiltonian can 



be found in [4]. Intuitively, since the first target time t\ 
is fixed, this choice generates the most direct traverse 
|^o) —> \(/>i), hence requiring least change in the Hamil- 
tonian at initial times t <C t\. 

The Euler-Lagrange equations governing stationary 
points of ([3| are obtained by taking the variation of J 
and requiring 5 J = 0. Writing A = ((5/7)£/ _1 we have 

5 J = I ^ ((iH, i6H) + (M, A — [UU-\A] + iSH) 

J t V 

1 m 

+(5M : UU- 1 + Ltf))dt + ^ E 5D2 W*i > 
= / ((M-ii^,i^) + (-M+[/7/7- 1 ,M],A) 
+ (<5M, f/C/" 1 + Lff>)d* + J2 SD 2 ^,^) 



j'=i 



m—i 

£ [(AMfo), A(t rf )> + (iAH( tj ),i5H( tj )) 



j = l 



+ (M(t m ),A(t m )) + (iH(t m ),i5H(t m )), 



(4) 



where in the second step we have integrated by parts, 
and used the notations AM(tj) = M(£~) - M(t+) and 
AJTfe) = JT(tT) - ff(tt), with M(£+) = lim Uti M(t) 
and M(tr) = lim nt . M(t); and similarly for H(tf). It 
follows from Q that on the open intervals j = 

0, . . . , m — 1, the following equations hold: 

Lff - M = 0, M + [M, = 0, + Lff = 0. (5) 

Additionally, at the nodes t = tj, we require matching 
conditions. To work them out, let us calculate the varia- 
tion 5D 2 = 2D5D appearing in Q. From the definition 
([2| and the relation 

^ (^|(1-£^)|0)(0|(1 + £^)|^) 



-eA 



(mm) , mm® (mm 



m)(4>\<t>) 



<</#)<V#> 



e + 0(e 2 ), (6) 



which holds for any A = -Al we find, bearing in mind 
that if D = 2arccos(v / a;) then AD /Ax = —2/ sin(£>), 



SD 



D(e EA x[>, 



Ae 



e=0 



-4jR[(V#)(0|A|V>)] 
Bin(U)<0|^)^|V> 



•(7) 



From (|7|), and writing ZX, = D($tji<t>j)i we deduce that 
^ 2 = 2Dj(Fj,A(tj)), where 



(8) 



The relevant matching conditions at the nodes are there- 
fore given by: 

ff(tt) - H(tT) = 0, M(tt) - M(iT) = DjFj/cr 2 , (9) 
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(a)cr = 0.04 



(b)cr = 0.01 



FIG. 2: (Colour online) Orbits on the state space generated by the 
solution to the quantum spline problem. The black dots indicate 
the initial (lower left) and the target points. The optimal trajecto- 
ries are shown for two different values of the tolerance parameter: 
a = 0.04 and a = 0.01. Lower values of the tolerance parameter 
translate, through the cost functional <J, into a stronger penalty 
on the mismatch. 



whereas we require H(t m ) = and M(t rn ) + D rn F rn / a 2 
at the terminal point. Quantum spline problems are 
therefore solved by finding a solution to equations ([5| and 
([9| that satisfies, in addition, the terminal conditions at 
t m . On open time intervals equation ([5| yields 



H + i[H,H] =0. 



(10) 



This is the right-reduced equation for the so-called Rie- 
mannian cubics on SU(n + 1) with respect to the bi- 
invariant Riemannian metric induced by the inner prod- 
uct Q. That is, U(t) is a Riemannian cubic on the open 
time intervals (^,^ + i). Here, by a Riemannian cubic we 
mean a solution to a certain fourth-order equation for 
a curve on a Riemannian manifold (see [7] for further 
details). The node conditions ([9| imply that U(t) is a 
Riemannian cubic spline, a twice continuously differen- 
tiable curve that is composed of a series of cubics. 

We remark on the important structure of the Lagrange 
multiplier M(t) implied by the equations of motion that 
makes it sufficient to consider a subspace of su(n + 1) 
when searching for the optimal initial value M(0). Let 
us denote by the totality of trace-free skew-Hermitian 
generators of unitary motions that leave the state in- 
variant, and Vj; its complement with respect to the in- 
ner product 0. Then, we have the following Lemma: 
M(t) G V^ t (this holds because for all j, DjFj G V^ t , 

and from (g), M(t) te{t Jtt = Ad um{tj+l) -i M(tJ +1 )). 
This result is significant, because the search for the 
optimal M(0) can be restricted to the 2n-dimensional 
subspace V^ Q of the n(n + 2)-dimensional Lie algebra 
5u(n + 1). 

Before we indicate the process for the implementation 
of the optimisation scheme, let us show some results first. 
Consider a two-level system (n = 1). We can think of this 
system as a spin-| particle immersed in a magnetic field. 
If n(t) is the unit direction of the field at time t, the 
Hamiltonian of the system can be written in the form 
H(t) = uo(t)cr -n(t), where uo(t) is the field strength. In 
this case the state space is just the Bloch sphere S 2 . We 





(a)Evolution of the rotation 
axis n(t) for a = 0.04 



(b) Evolution of the rotation 
axis n(t) for a = 0.01 
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(c)Field strength uo(t) for 
a = 0.04 
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(d)Field strength uj(t) for 
a = 0.01 



FIG. 3: (Colour online) The quantum spline H(t). Hamiltonians 
that generate the dynamical trajectories in Fig. [2] The top row 
shows the orbits of the endpoint of the rotation axis n(t). The 
bottom row shows the field strength uo(t). These images illustrate 
the fact that as the value of a is decreased, the amount of change 
in the optimal Hamiltonian Hit) increases. 



have implemented the optimisation for a set of target 
states on £ 2 , an initial state |^o)j an d a set of times. 
Using the resulting Hamiltonian we have generated the 
dynamics of the state, as illustrated in Fig. [l] In Fig. [2] 
we have sketched the effect of choosing different tolerance 
levels. When the value of a is reduced, the resulting orbit 
\?pt} traverses closer to the vicinities of the target states 
{|^)}. FromQ, one sees that this may be realised at the 
expense of varying the Hamiltonian H(t) more rapidly. 
This effect can be visualised in the case of a two-level 
system, since H{t) is characterised by the the scalar field 
strength uo(t) and the unit vector n(t) G R 3 . In Fig. [3] 
we have plotted the end-point of the unit vector n(t) on 
a sphere, and the values of cj(t), for different choices of 
a. These plots show that both n(t) and uo(t) vary more 
rapidly at smaller tolerance level (i.e. smaller a). 

Another example we consider here is a controlled mo- 
tion of a quantum state on the coherent-state subspace of 
the state space. Consider SU(n + l) coherent states [HJ[9] 
in arbitrary dimensions. These coherent states can be 
generated by taking symmetric tensor products of 'single- 
particle' states. In the context of quantum information 
theory, these states correspond to totally disentangled 
states inside the symmetric subspace of the Hilbert space 
of the combined system. Each coherent state thus cor- 
responds to the image of a map, known as the Veronese 
embedding [TQl [TT], of a pure state. Therefore, given a 
set of points on a coherent-state space we identify them 
with states on a single-particle Hilbert space, solve the 
quantum spline problem as indicated above, and map the 
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result back to the larger Hilbert space. In particular, the 
coherent quantum spline is generated by the symmet- 
ric tensor product Hamiltonian ® s H(t). This elemen- 
tary procedure works because (a) the Veronese embed- 
ding commutes with the action of SU(n + 1); and (b) 
the natural metrics on the spaces of coherent states are 
scalar-multiples of the metric Q used here [TO] . 

Next we discuss a numerical approach for finding a lo- 
cal minimum of the cost functional (p3). The search can 
be restricted to solutions of (|5| and (9|, which are en- 
coded by their initial conditions M(to) and H(to). We 
can therefore regard ([3| as a function of these initial con- 
ditions, and perform a descent algorithm on that func- 
tion. The terminal conditions at t m can then be used to 
test whether we have arrived at a local minimum. 

For a numerical implementation we can discretise the 
equations of motion ([5| and ([9|, and find the approximate 
gradient of J\ alternatively, we can introduce an approxi- 
mation Jd of J defined on a discrete path space, and take 
its variation, which yields a set of discrete equations of 
motion. Here we follow the latter method, which permits 
the use of adjoint equations [6 for an efficient calculation 
of the exact gradient of Jd- This method is highly effec- 
tive in dealing with higher-dimensional (n > 1) systems. 
Moreover, in this method discrete critical curves of Jd 
satisfy a version of the terminal conditions at t = t m 
exactly, and this leads to a precise method for testing 
convergence. In addition, such curves fulfil the condi- 
tions for the above-stated Lemma on their discrete time 
domain, which can be exploited by restricting the search 
for the optimal initial value of M to V^ Q . 

The implementation will make use of the Cayley map 
r : 5u(n + 1) — >• SU(n + 1), which approximates the 
Lie exponential according to I 4 (1 - X/2) _1 (l + 
X/2). We will also need the left-trivialised differential d\\ 
dir x Y = (d/de)r(X + 6:y')r(X)" 1 | £= o, which is given by 
(1 - X/2)~ 1 Y(t + X/2)- 1 . We discretise the time in- 
terval t m — to into N steps such that (t m — to)/N = h, 
and we let t^ = to + fih for fi = 0, . . . , N. For simplic- 
ity, we assume that the nodal times {t,}j=o,...,m coin- 
cide with some of the discrete time steps t nj = to + rijh, 
where no = and n m = N. To obtain a discrete ver- 
sion of the cost functional, we approximate the time 
derivative — \H of the generator by the discrete vari- 
ables {Lj^}. The complete set of discrete variables is 
therefore (U^, ii^, M M , with ji = 0, . . . , N. Writing 
A M = 5^ n .DjFj/a 2 and making use of the Euler method 
of [6] , we obtain the following set of discrete equations of 
motion for /i = 0, . . . , N — 1: 



M M+ i = K^ +i )(d i r_ i ^ +1 )(M M + A /i ) 



= L, 



h(din hHll+1 )M^ +1 



(11) 



U, 



= ri-ihH^U, 



H^+i = Hp 



ihL 



These equations can be integrated for given initial values 
Mo and Lq (recall that Uq = 1 and Ho are prescribed). 



The terminal conditions are Ln = and M/v + A at = 
0. The discrete cost functional Jd in terms of initial 
conditions (Mq,Lq) is 



N-l 



h 



1 

2^ 



^Z) 2 (/7 n >, (12) 



whereby the equations of motion (11) are implied. 



A local minimum can be found by a gradient descent 
method, which requires the computation of the gradient 
of Jd- The estimation of the gradient via finite- difference 
methods requires the repeated forward integration of the 
system of equations (11). The number of forward in- 
tegrations increases with the number of dimensions of 
the Lie algebra (n 2 to leading order). Such estimation 
procedures thus quickly become unfeasible for higher- 
dimensional systems. This difficulty can be avoided by 
using the method of adjoint equations, which can be 



readily implemented for the discretisation (11), (12) pre- 



sented here (see Supplemental Material for details, and 
arXiv: 1206. 2675 v2 for a numerical code). Then, the ex- 
act gradient is obtained at the cost of integrating twice 
(once forward, once backward) a system of equations of 
the same complexity as (11). This allows for an efficient 



treatment of the quantum spline problem when n > 1. 
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Gradient computation via adjoint equations: 
Supplemental Material 

Here we describe the method of adjoint equations for 
the efficient computation of the gradient of We will 
supply the necessary equations, referring to [6] for fur- 
ther details. First we need convenient expressions for 
the partial derivatives of Jd with respect to the variables 
Mo and L$. These are denoted by V M Q Jd and V ' L Q Jd, 
respectively, and are defined by 



SJd 



de 



J d (M + eSM ,Lo + eSL ) 



e=0 



(V Mo Jd,5M ) + (V Lo J d ,5L ) , 



(13) 



for all 5M , SLq in su(n + 1). 

Besides the Cayley map r and the left-trivialised differ- 
ential d\ we will also need the right-trivialised differential 
d r : d r rxY = T(X)~ 1 (d/de)r(X-\-sY)\ £= o, which is given 
by (1 + X/2)- 1 F(l - X/2)~\ We define the functions 
Fj for j = 1, . . . m: 



sin(I)(^0)<^><V#> ■ 



as well as the functions A M for \i = 0, . . . , TV defined by 
A fl ^) = S^ 3 D(^,4> j )F j ^)/a 2 . 

Next, define an augmented functional Q, in which the 
discrete equations of motion ( pTj ) are incorporated us- 
ing Lagrange multipliers. These Lagrange multipliers 
will be denoted {P^P^V^V^}^...^. Let us intro- 
duce the shorthand notation x representing the discrete 
path {C/^, iii^, M^, £ / x} / z=o,... J iv an d ^ representing the 
Lagrange multipliers {P^, V^, V^f} M= i 5 ... ? iv- Writing 
IVv) = k^l^o) the augmented functional is given by 



AT-l 



0(*,A) 



i(L / ,,^) + (P M +1 ,r- 1 (^ + iC/- 1 ) 

/x=0 



+(^ +1 ,Wt Mm+1 )M m+1 



1 m 



2a : 



No constraints are assumed here, apart from the pre- 
scribed initial Hamiltonian Ho and £7(0) = 1. Note that 
Q(x : A) = Jd(Mo,Lo) for any choice of Lagrange multi- 
pliers A, provided that x satisfies the discrete equations 



of motion (11) for given initial values Mq and Lq . More- 



over, taking variations of Q yields 
SO = h(L -hP}-V?,5L ) - h(d r T- ihHl V?,5Mo) (14) 



if x satisfies (11) and A is a solution of the adjoint equa- 
tions. To specify the adjoint equations, we introduce 



functions K ± by the defining relation 
d 



de 



(dir ±hix +eY)M, V) 



£ = 



= (Kfx,M) V > Y ) 



for all X, Y, M, V in su(n + 1), and define by 



-(A, (e"V),V) 



M M (V-,v),A> 



e=0 



for all V,ie su(n + l) and state vectors The adjoint 
equations consist of conditions at the final time point: 



V° 
v N 

pO 



0, 



0, 



-(diT- ihHN )A N (ip N ), 



(15) 
(16) 



and the following set of equations 
V° = V° +1 +hP' +1 -L, 

K = d r T ihH„ [drT-ihH^V^ - hd r T ihHtl V°] 

dirZl hH ^ + Pl +1 ~ \A^) 

+ A^(lp M d r T- ihHll+1 V^ +1 ) 



P° = dm hHiJ 



(17) 



P 1 — P 1 _l h p0 _ h j{ yO _ is- V 1 

for /i = 1, . . . , TV — 1. These equations are posed back- 
wards. That is, solving the adjoint equations entails ini- 
tialising the Lagrange multipliers at time point TV ac- 



cording to (15) and (16), and then iterating backwards 



from \i = TV to /i = 1 using (17). 



We now obtain the exact gradient of J& from (14). In- 
deed, let (Mo(e), Lo(e)) be variations of initial conditions 
(Mo,i>o), let x(e) be the corresponding solutions to the 
discrete equations of motion (11) with x = x(0), and let 
A be a solution to the adjoint equations (Il]^-([l7|). Then 



SJ d = ^Jd(Mo(e),L (s)) 
de 



£ = 



de 



£ = 



h(L - hPl - V^SLo) - h(drT- ihHl V?,5Mo), 



where we have used (14) in the last equality. Recalling 



definition (|13|) we conclude that 



hDr- ihHl V?, V Lo J d = h(L - hP\ - V?) . 

(18) 



In summary, the partial derivatives V M Jd and V L Q Jd 
are computed as follows: (i) Integrate the system of equa- 
tions (11 ) up to /i = TV; (ii) Initialise the Lagrange multi- 
pliers at \i = TV according to (15) and (16); (hi) Integrate 
backwards the system of equations (17) until /i = 1; and 



(iv) Obtain V M Q Jd and V L Jd by evaluating the right 
hand sides of (18). 



